###
# Replication code for 
# Rallying around the flag in times of COVID-19
###

require(tidyverse)
require(magrittr)
require(rdd)
require(here)
require(fastDummies)
require(knitr)

#### Data ####

# NB: Make sure the data file is stored in your active project's root folder. 
# Otherwise, change the use of here::here.

data <- read.csv(here::here("replication_data_covid_19.csv")) %>% as.tibble
names(data) <- gsub("\\."," ",names(data))

###
# Manuscript
###

#### Figure 1 ####

daily_means <- data %>%
  rename(`Trust in PM's\nAdministration` = `Trust in Government`) %>%
  mutate(time = round(numeric_time,0) - .5) %>%
  select(`Trust in PM's\nAdministration`,time) %>%
  group_by(time) %>%
  summarize(
    daily_mean = mean(`Trust in PM's\nAdministration`,na.rm = T)
  ) %>% rename(numeric_time = time)

data %>%
  select(contains("Trust in Government"),numeric_time) %>%
  gather(variable,value,-numeric_time) %>%
  ggplot(aes(x = numeric_time, y = value, group = numeric_time >= 0)) +
  geom_point(color = "lightgray",alpha = .4,size = .7) + 
  geom_smooth(formula = y ~ x,se=T,color="black",method = "lm",alpha = .3) +
  theme_minimal() +
  theme(
    panel.spacing = unit(2,"lines"),
    strip.text.y = element_text(size = 8, angle = 0)
  ) + 
  scale_y_continuous(limits=c(0,10),expand = c(0,0)) +
  scale_x_continuous(breaks = seq(-2,20,by=2))+ 
  geom_vline(xintercept = 0) + 
  labs(
    x = "Time from lockdown\n(in days)",
    y = "Trust in PM's administration"
  ) + 
  geom_point(data=daily_means,aes(x = numeric_time, y = daily_mean))

#### Table 1 ####

model1 <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0))
model1_tidy <- model1 %>% broom::tidy() %>% mutate(model = 1)

model2 <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age)
model2_tidy <- model2 %>% broom::tidy() %>% mutate(model = 2)

model3 <- lm(data=data,
             `Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
               Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
               Children)
model3_tidy <- model3 %>% broom::tidy() %>% mutate(model = 3)

model4 <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
               Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
               Children + Expectation)
model4_tidy <- model4 %>% broom::tidy() %>% mutate(model = 4)

models <- bind_rows(model1_tidy,model2_tidy,model3_tidy,model4_tidy) %>%
  group_by(model) %>%
  mutate_if(is.numeric,formatC,digits=3,format="f") %>%
  mutate(id = row_number()) %>%
  mutate(estimate = paste0(estimate," (",std.error,")")) %>%
  select(id,term,estimate,model) %>%
  spread(model,estimate) %>%
  select(-id) %>%
  mutate(term = gsub("\\`","",ifelse(grepl(">=",term),"Lockdown",term))) %>%
  transmute_all(funs(. = ifelse(is.na(.),"",.)))

names(models) <- c("","model 1","model 2","model 3","model 4")

models

#### Figure 2 ####

figure2_data <- data %>%
  rename(`Trust in PM's\nAdministration` = `Trust in Government`) %>%
  gather(variable,value,contains("Trust")) %>%
  mutate(
    trust = value,
    variable = variable,
    numeric_time = numeric_time
  ) %>%
  group_by(variable) %>%
  nest() %>%
  mutate(rdd = map(.x = data,~rdd::RDestimate(data=.x,trust~numeric_time,cutpoint = 0,bw = 7)))

figure2_its <- figure2_data %>%
  mutate(est = map_dbl(rdd,~.x$est[["LATE"]])) %>%
  mutate(se = map_dbl(rdd,~.x$se[[1]])) %>%
  mutate(p_value = 2*pnorm(-abs(est/se))) %>%
  mutate(
    lower = est - 1.96*se,
    upper = est + 1.96*se,
    Lower = est - 1.6*se,
    Upper = est + 1.6*se
  ) %>% 
  mutate(estimator = "ITS") %>%
  select(-data,-rdd) %>% ungroup()

figure2_ols <- figure2_data %>%
  mutate(rdd = map(data,~lm(data=.,trust~I(numeric_time >= 0)) %>% broom::tidy(.))) %>%
  unnest(rdd) %>%
  filter(term == "I(numeric_time >= 0)TRUE") %>%
  mutate(
    se = std.error,
    est = estimate,
    lower = est - 1.96*se,
    upper = est + 1.96*se,
    Lower = est - 1.6*se,
    Upper = est + 1.6*se,
    p_value = p.value,
    estimator = "OLS"
  ) %>%
  select(variable,est,se,p_value,lower,upper,Lower,Upper,estimator) %>% ungroup()

figure2_ols_controls <- figure2_data %>%
  mutate(rdd = map(data,~lm(data=.,trust~I(numeric_time >= 0) + Gender + Age + Education + `Employment status` + 
                              `Unemployment duration` + `Unemployment periods` + Children + 
                              `Too many activities` + `Too many demands` + `Too little control`) %>% 
                     broom::tidy(.))) %>%
  unnest(rdd) %>%
  filter(term == "I(numeric_time >= 0)TRUE") %>%
  mutate(
    se = std.error,
    est = estimate,
    lower = est - 1.96*se,
    upper = est + 1.96*se,
    Lower = est - 1.6*se,
    Upper = est + 1.6*se,
    p_value = p.value,
    estimator = "OLS w. controls"
  ) %>%
  select(variable,est,se,p_value,lower,upper,Lower,Upper,estimator) %>% ungroup()

plot_coefficients <- rbind(
  figure2_ols,
  figure2_ols_controls,
  figure2_its
) %>% mutate(variable = factor(variable, levels = rev(c("Trust in PM's\nAdministration",
                                                        "Trust in Parliament",
                                                        "Trust in the Public Sector",
                                                        "Trust in Courts",
                                                        "Trust in the Media")))) %>%
  mutate(estimator = factor(estimator, levels = c("OLS","OLS w. controls","ITS")))


coefficient_plot <- plot_coefficients %>%
  arrange(variable) %>%
  ggplot(aes(x = variable, y = est, ymin = lower, ymax = upper, shape = estimator)) +
  geom_point(position = position_dodge(.5),size=3) +
  geom_errorbar(width = 0,position = position_dodge(.5)) +
  geom_errorbar(data=plot_coefficients,aes(ymin = Lower, ymax = Upper),width=0,position = position_dodge(.5),lwd=1.2) + 
  geom_hline(yintercept = 0, lty = 4, alpha = .5) +
  scale_y_continuous(expand = c(.05,0,0,.2),breaks = seq(-.5,2,length.out = 6)) + 
  theme_minimal() + 
  theme(
    legend.position = "bottom",
    axis.text = element_text(size = 12),
    legend.text = element_text(size = 12)
  ) +
  labs(
    x = "Dependent variable",
    y = "Estimate",
    shape = "Estimator: "
  ) + coord_flip()

coefficient_plot + theme(
  axis.text.x = element_text(size = 16),
  axis.text.y = element_text(size = 16),
  axis.title.x = element_text(size = 18),
  axis.title.y = element_text(size = 18),
  legend.text = element_text(size = 16)
)

#### Table 2 ####

figure2_its %>% 
  mutate(`Dependent variable` = ifelse(str_detect(variable,"Administration"),"Trust in PM's Administration",variable)) %>%
  select(`Dependent variable`,est,se,p_value) %>%
  mutate_if(is.numeric,formatC,digits = 3, format = "f") %>%
  knitr::kable()

###
# Supplementary Information
###

#### SI section 3: Descriptive statistics ####

data %>%
  select(`Trust in Government`:`Too little control`) %>%
  rename(`Trust in PM's Administration` = `Trust in Government`) %>%
  fastDummies::dummy_cols(select_columns = "Gender",remove_first_dummy = F) %>%
  fastDummies::dummy_cols(select_columns = "Employment status",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Employment status")),funs(ifelse(`Employment status_NA` == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Unemployment duration",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Unemployment duration")),funs(ifelse(`Unemployment duration_NA` == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Children",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Children")),funs(ifelse(Children_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Education",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Education")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  select(-Gender,-`Employment status`,-`Unemployment duration`,-Children,-Education,-contains("_NA")) %>%
  gather(variable,value) %>%
  group_by(variable) %>%
  summarize(
    N = sum(!is.na(value)),
    `Item response rate` = N / length(value),
    Mean = mean(value,na.rm = T),
    `Std.dev` = sd(value,na.rm = T),
    Min = min(value,na.rm = T),
    Max = max(value,na.rm = T)
  ) %>%
  arrange(desc(str_detect(variable,"Trust"))) %>%
  mutate(variable = str_replace_all(variable,"_",": ")) %>%
  mutate_at(vars(`Item response rate`,Mean,`Std.dev`),funs(formatC(.,digits = 3,format = "f"))) %>%
  knitr::kable(.)

#### SI section 4: Developments over time #### 

data %>%
  select(contains("Trust"),numeric_time) %>%
  gather(variable,value,-numeric_time) %>%
  ggplot(aes(x = numeric_time, y = value, group = numeric_time >= 0)) +
  geom_point(color = "lightgray",alpha = .4,size = .7) + 
  geom_smooth(formula = y ~ x,se=T,color="black",method = "lm") +
  facet_grid(variable~.) + 
  theme_minimal() +
  theme(
    panel.spacing = unit(2,"lines"),
    strip.text.y = element_text(size = 8, angle = 0)
  ) + 
  scale_y_continuous(limits=c(0,10),expand = c(0,0)) +
  scale_x_continuous(breaks = seq(-2,20,by=2))+ 
  geom_vline(xintercept = 0) + 
  labs(
    x = "Time from lockdown\n(in days)",
    y = "Trust"
  )

#### SI section 5: Balance tests ####

gender_p <- chisq.test(with(data,table(Gender,numeric_time >= 0)))$p.value
age_p <- broom::tidy(lm(data=data,Age ~ I(numeric_time >= 0)))$p.value[[2]]
education_p <- chisq.test(with(data,table(Education,numeric_time >= 0)))$p.value
status_p <- chisq.test(with(data,table(`Employment status`,numeric_time >= 0)))$p.value
duration_p <- chisq.test(with(data,table(`Unemployment duration`,numeric_time >= 0)))$p.value
periods_p <- broom::tidy(lm(data=data,`Unemployment periods` ~ I(numeric_time >= 0)))$p.value[[2]]
children_p <- chisq.test(with(data,table(`Children`,numeric_time >= 0)))$p.value

data %>%
  fastDummies::dummy_cols(select_columns = "Gender",remove_first_dummy = F) %>%
  fastDummies::dummy_cols(select_columns = "Education",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Education")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Children",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Children")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Employment status",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Employment status")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Unemployment duration",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Unemployment duration")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  select(contains("Gender_"),contains("Education_"),
         contains("Employment status_"),contains("Unemployment duration_"),
         contains("Children_"),
         Age,`Unemployment periods`,numeric_time) %>%
  mutate(lockdown = numeric_time >= 0) %>%
  select(-contains("_NA"),-numeric_time) %>%
  group_by(lockdown) %>%
  summarize_all(
    funs(mean(.,na.rm=T))
  ) %>% na.omit %>%
  gather(variable,value,-lockdown) %>%
  mutate(lockdown = ifelse(lockdown == T,"After lockdown","Before lockdown")) %>%
  group_by(variable) %>%
  spread(lockdown,value) %>%
  ungroup() %>%
  mutate_at(vars(contains("lockdown")),funs(round(.,3))) %>%
  mutate(variable = str_replace_all(variable,"_",": ")) %>%
  mutate(Order = c(1,5,2,3,4,7,6,9,8,11,10,12,13,14,16,17,18,15,19)) %>%
  arrange(Order) %>% select(-Order) %>%
  mutate(`P value` = c(round(age_p,3),
                       round(children_p,3),"","","",
                       round(education_p,3),"","","",
                       round(status_p,3),"",
                       round(gender_p,3),"",
                       formatC(duration_p,digits = 3,format = "f"),"","","","",
                       formatC(periods_p,digits = 3, format = "f"))) %>% 
  knitr::kable()

reduced_data <- data %>% filter(abs(numeric_time) <= 7)

gender_p <- chisq.test(with(reduced_data,table(Gender,numeric_time >= 0)))$p.value
age_p <- broom::tidy(lm(data=reduced_data,Age ~ I(numeric_time >= 0)))$p.value[[2]]
education_p <- chisq.test(with(reduced_data,table(Education,numeric_time >= 0)))$p.value
status_p <- chisq.test(with(reduced_data,table(`Employment status`,numeric_time >= 0)))$p.value
duration_p <- chisq.test(with(reduced_data,table(`Unemployment duration`,numeric_time >= 0)))$p.value
periods_p <- broom::tidy(lm(data=reduced_data,`Unemployment periods` ~ I(numeric_time >= 0)))$p.value[[2]]
children_p <- chisq.test(with(reduced_data,table(`Children`,numeric_time >= 0)))$p.value

reduced_data %>%
  fastDummies::dummy_cols(select_columns = "Gender",remove_first_dummy = F) %>%
  fastDummies::dummy_cols(select_columns = "Education",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Education")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Children",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Children")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Employment status",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Employment status")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  fastDummies::dummy_cols(select_columns = "Unemployment duration",remove_first_dummy = F) %>%
  mutate_at(vars(contains("Unemployment duration")),funs(ifelse(Education_NA == 1,NA,.))) %>%
  select(contains("Gender_"),contains("Education_"),
         contains("Employment status_"),contains("Unemployment duration_"),
         contains("Children_"),
         Age,`Unemployment periods`,numeric_time) %>%
  mutate(lockdown = numeric_time >= 0) %>%
  select(-contains("_NA"),-numeric_time) %>%
  group_by(lockdown) %>%
  summarize_all(
    funs(mean(.,na.rm=T))
  ) %>% na.omit %>%
  gather(variable,value,-lockdown) %>%
  mutate(lockdown = ifelse(lockdown == T,"After lockdown","Before lockdown")) %>%
  group_by(variable) %>%
  spread(lockdown,value) %>%
  ungroup() %>%
  mutate_at(vars(contains("lockdown")),funs(round(.,3))) %>%
  mutate(variable = str_replace_all(variable,"_",": ")) %>%
  mutate(Order = c(1,5,2,3,4,7,6,9,8,11,10,12,13,14,16,17,18,15,19)) %>%
  arrange(Order) %>% select(-Order) %>%
  mutate(`P value` = c(round(age_p,3),
                       formatC(children_p,digits = 3,format = "f"),"","","",
                       round(education_p,3),"","","",
                       round(status_p,3),"",
                       round(gender_p,3),"",
                       formatC(duration_p,digits = 3,format = "f"),"","","","",
                       formatC(periods_p,digits = 3, format = "f"))) %>% 
  knitr::kable()

#### SI section 6: OLS restricted samples ####

data %>%
  nest() %>%
  crossing(include_up_to = seq(2,20,length.out = 100)) %>%
  mutate(restricted_data = map2(.x = data, .y = include_up_to,~filter(.x,abs(numeric_time) < .y))) %>%
  mutate(model = map(restricted_data,~lm(data=.,
                                         `Trust in Government` ~ I(numeric_time >= 0) + 
                                           Gender + Age + Education + `Employment status` +
                                           `Unemployment duration` + `Unemployment periods` +
                                           Children + Expectation))
  ) %>%
  mutate(tidy_model = map(model,~broom::tidy(.,conf.int = T))) %>%
  unnest(tidy_model) %>%
  filter(term == "I(numeric_time >= 0)TRUE") %>%
  ggplot(aes(x = include_up_to, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_line() + 
  geom_point(color = "lightgray",alpha = .4,size = .7) + 
  geom_ribbon(alpha = .3) + 
  theme_minimal() +
  theme(
    panel.spacing = unit(2,"lines"),
    strip.text.y = element_text(size = 8, angle = 0)
  ) + 
  scale_y_continuous(limits=c(-.5,1.75)) +
  scale_x_continuous(breaks = seq(2,20,by=2),expand = c(0,0,0,.4)) + 
  geom_hline(yintercept = 0,lty=3) + 
  labs(
    x = "Data included\n(in days from lockdown)",
    y = "OLS estimate for lockdown"
  )

#### SI section 7: ITS bandwidth sensitivity ####

its_models <- data %>%
  rename(`Trust in PM's\nAdministration` = `Trust in Government`) %>%
  gather(variable,value,contains("Trust")) %>%
  mutate(
    trust = value,
    variable = variable,
    numeric_time = numeric_time
  ) %>%
  group_by(variable) %>%
  nest() %>%
  crossing(bw = seq(1,20,length.out = 100)) %>%
  mutate(rdd = map2(.x = data, .y = bw,~rdd::RDestimate(data=.x,trust~numeric_time,cutpoint = 0,bw = .y)))

its_models %>%
  mutate(est = map_dbl(rdd,~.x$est[["LATE"]])) %>%
  mutate(se = map_dbl(rdd,~.x$se[[1]])) %>%
  mutate(
    lower = est - 1.96*se,
    upper = est + 1.96*se
  ) %>% 
  mutate(variable = ifelse(variable == "Trust in the Public Sector","Trust in the\nPublic Sector",variable)) %>%
  ggplot(aes(x = bw, y = est, ymin = lower, ymax = upper, group = variable)) +
  geom_hline(yintercept = 0,alpha = .4,lty = 4) +
  geom_ribbon(alpha = .3) + 
  geom_line() + 
  facet_grid(.~variable) + 
  scale_x_continuous(breaks = seq(0,20,by = 5),limits = c(0,20)) +
  theme_minimal() + 
  theme(axis.text.x = element_text(size = 8),
        panel.spacing = unit(2, "lines"),
        plot.margin = margin(0, 1, 1, 1, "cm"),
        strip.text = element_text(size = 8, face = "bold",angle = 90)) + 
  labs(
    x = "Discontinuity Bandwidth",
    y = "LATE Estimate"
  )

#### SI section 8: ADCE ####

first_stage_activities <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
   Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
   Children + Expectation + `Too many activities`)

direct_activities <- lm(data=data,
  I(`Trust in Government` - `Too many activities`*coef(first_stage_activities)["`Too many activities`"]) ~
    I(numeric_time >= 0) + Gender + Age + 
    Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
    Children + Expectation
  )

first_stage_demands <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
   Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
   Children + Expectation + `Too many demands`)

direct_demands <- lm(data=data,
    I(`Trust in Government` - `Too many demands`*coef(first_stage_demands)["`Too many demands`"]) ~
      I(numeric_time >= 0) + Gender + Age + 
      Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
      Children + Expectation
)

first_stage_control <- lm(data=data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
                            Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                            Children + Expectation + `Too little control`)

direct_control <- lm(data=data,
                     I(`Trust in Government` - `Too little control`*coef(first_stage_control)["`Too little control`"]) ~
                       I(numeric_time >= 0) + Gender + Age + 
                       Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                       Children + Expectation
)

adce <- c(
  coef(direct_activities)[2],
  coef(direct_demands)[2],
  coef(direct_control)[2]
)

boots <- matrix(nrow = 2000,ncol = 3)

for(i in 1:nrow(boots)){
  sample_data <- data[sample(1:nrow(data),replace = TRUE),]
  
  first_stage_activities <- lm(data=sample_data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
                                 Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                                 Children + Expectation + `Too many activities`)
  
  direct_activities <- lm(data=sample_data,
                          I(`Trust in Government` - `Too many activities`*coef(first_stage_activities)["`Too many activities`"]) ~
                            I(numeric_time >= 0) + Gender + Age + 
                            Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                            Children + Expectation
  )
  
  first_stage_demands <- lm(data=sample_data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
                              Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                              Children + Expectation + `Too many demands`)
  
  direct_demands <- lm(data=sample_data,
                       I(`Trust in Government` - `Too many demands`*coef(first_stage_demands)["`Too many demands`"]) ~
                         I(numeric_time >= 0) + Gender + Age + 
                         Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                         Children + Expectation
  )
  
  first_stage_control <- lm(data=sample_data,`Trust in Government` ~ I(numeric_time >= 0) + Gender + Age + 
                              Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                              Children + Expectation + `Too little control`)
  
  direct_control <- lm(data=sample_data,
                        I(`Trust in Government` - `Too little control`*coef(first_stage_control)["`Too little control`"]) ~
                          I(numeric_time >= 0) + Gender + Age + 
                          Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                          Children + Expectation
  )
  
  boots[i,] <- c(
    coef(direct_activities)[2],
    coef(direct_demands)[2],
    coef(direct_control)[2]
  )
}

adce.se <- apply(boots,2,"sd")

tibble(
  `Blocked Mediator variable` = c("Too many activities","Too many demands","Too little control"),
  `ADCE Lockdown Estimate` = adce,
  `SE` = adce.se
) %>%
  mutate(
    statistic = `ADCE Lockdown Estimate` / SE,
    `P value` = 2*pnorm(-abs(statistic))
  ) %>% select(-statistic) %>%
  mutate_if(is.numeric,formatC,digits = 3, format = "f") %>%
  mutate(`ADCE Lockdown Estimate` = paste0(`ADCE Lockdown Estimate`," (",SE,")")) %>%
  select(-SE) %>%
  knitr::kable()
  
#### SI section 9 Moving discontinuity ####

its_models_threshold <- data %>%
  rename(`Trust in PM's\nAdministration` = `Trust in Government`) %>%
  gather(variable,value,contains("Trust")) %>%
  mutate(
    trust = value,
    variable = variable,
    numeric_time = numeric_time
  ) %>%
  group_by(variable) %>%
  nest() %>%
  crossing(cutoff = seq(-1,5,length.out = 100)) %>%
  mutate(rdd = map2(.x = data, .y = cutoff,~rdd::RDestimate(data=.x,trust~numeric_time,cutpoint = .y,bw = 7)))

its_models_threshold %>%
  mutate(est = map_dbl(rdd,~.x$est[["LATE"]])) %>%
  mutate(se = map_dbl(rdd,~.x$se[[1]])) %>%
  mutate(
    lower = est - 1.96*se,
    upper = est + 1.96*se
  ) %>% 
  mutate(variable = ifelse(variable == "Trust in the Public Sector","Trust in the\nPublic Sector",variable)) %>%
  ggplot(aes(x = cutoff, y = est, ymin = lower, ymax = upper, group = variable)) +
  geom_hline(yintercept = 0,alpha = .4,lty = 4) +
  geom_ribbon(alpha = .3) + 
  geom_line() + 
  facet_grid(.~variable) + 
  geom_vline(xintercept = 0) +
  scale_x_continuous(breaks = seq(-1,5,by = 1)) +
  theme_minimal() + 
  theme(axis.text.x = element_text(size = 8),
        panel.spacing = unit(2, "lines"),
        plot.margin = margin(0, 1, 1, 1, "cm"),
        strip.text = element_text(size = 8, face = "bold",angle = 90)) + 
  labs(
    x = "Change in Discontinuity from lockdown\n(in days)",
    y = "LATE Estimate"
  )

ols_models_threshold <- data %>%
  rename(`Trust in PM's\nAdministration` = `Trust in Government`) %>%
  gather(variable,value,contains("Trust")) %>%
  mutate(
    trust = value,
    variable = variable,
    numeric_time = numeric_time
  ) %>%
  group_by(variable) %>%
  nest() %>%
  crossing(cutoff = seq(-1,5,length.out = 100)) %>%
  mutate(ols = map2(.x = data, .y = cutoff,~lm(data=.x,trust~I(numeric_time >= .y) + Gender + Age + 
                                                 Education + `Employment status` + `Unemployment duration` + `Unemployment periods` +
                                                 Children + Expectation)))

ols_models_threshold %>%
  mutate(tidy_ols = map(ols,~broom::tidy(.) %>% filter(str_detect(term,"numeric_time")))) %>%
  unnest(tidy_ols) %>%
  mutate(
    lower = estimate - 1.96*std.error,
    upper = estimate + 1.96*std.error
  ) %>% 
  mutate(variable = ifelse(variable == "Trust in the Public Sector","Trust in the\nPublic Sector",variable)) %>%
  ggplot(aes(x = cutoff, y = estimate, ymin = lower, ymax = upper, group = variable)) +
  geom_hline(yintercept = 0,alpha = .4,lty = 4) +
  geom_ribbon(alpha = .3) + 
  geom_line() + 
  facet_grid(.~variable) + 
  geom_vline(xintercept = 0) +
  scale_x_continuous(breaks = seq(-1,5,by = 1),limits = c(-1,5)) +
  theme_minimal() + 
  theme(axis.text.x = element_text(size = 8),
        panel.spacing = unit(2, "lines"),
        plot.margin = margin(0, 1, 1, 1, "cm"),
        strip.text = element_text(size = 8, face = "bold",angle = 90)) + 
  labs(
    x = "Change in Discontinuity from lockdown\n(in days)",
    y = "OLS Estimate"
  )

